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The XY model with quenched random phase shifts is studied by a T = finite size defect energy 
scaling method in 2d and 3d. The defect energy is defined by a change in the boundary conditions 
from those compatible with the true ground state configuration for a given realization of disorder. 
A numerical technique, which is exact in principle, is used to evaluate this energy and to estimate 
the stiffness exponent 9. This method gives 9 = —0.36 ± 0.013 in 2d and 9 = +0.31 ± 0.015 in Sd, 
which are considerably larger than previous estimates, strongly suggesting that the lower critical 
dimension is less than three. Some arguments in favor of these new estimates are given. 
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The XY model with random quenched disorder, also 
known as a gauge glass, as a model for a superconducting 
glass phase has been intensively investigated over the last 
decade. Since the mean field theory flux lattice moves in 
response to a transport current, the existence of disorder, 
which, however, destroys the flux lattice structure, be- 
comes essential in order to have a superconducting phase 
in a high superconductor. Since the pinning of flux 
lines by random disorder can prevent dissipation, a su- 
percurrent can exist in the glass phase 

From numerical [^-^ and experimental |^ studies, it 
is believed that the gauge glass has no ordered phase at 
any finite temperature in two dimensions. In three di- 
mensions, numerical studies indicate that the lower crit- 
ical dimension seems to be close to three. However the 
situation is less conclusive, since the simulations are lim- 
ited to small system sizes. Experimentally there is also 
some evidence which supports a finite temperature phase 
transition to a superconducting glass phase . 

The Hamiltonian of the gauge glass model is given by 

H=Y^V{e.,-6,-M,) (1) 

where V{4i) is a 27r periodic function of with a 
maximum at = tt, usually taken to be = 
—JijCos{4>ij). 6i is the phase of the superconducting or- 
der parameter at site i of a square lattice in 2d and a sim- 
ple cubic lattice in Zd. The sum is taken over all nearest 
neighbor sites. The random bond variables, Aij, which 
are responsible for the random frustration, are uniformly 
distributed in (— tt, tt]. The coupling constants, Jy, are 
taken as constant, — J. In this paper, we use a 
domain wall renormalization group (DWRG) ||l^,|l^, or 
defect energy scaling, method to investigate the possibil- 
ity of an ordered phase at small but finite temperature 
T. The idea is to find the ground state (GS) energy 
Eq{L) of a system of linear size L, then to change the 
boundary conditions (BC) in one direction to introduce 
a defect or domain wall into the system and finally to 
find the energy Ed{L) of this. The energy difference 



AE{L) = Ed{L) — Eq{L) is then the energy of a domain 
wall (defect) of length L. This is assumed to scale as 

< A£;(L) >- (2) 

where < ■ • • > means an average over disorder and 9 
is the stiffness exponent. The sign of 9 determines if 
a low T ordered phase exists or not. If < 0, in the 
thermodynamic limit the defect energy vanishes and at 
any T > will proliferate and destroy the order. On the 
other hand, if 6* > 0, the defect energy diverges as L — > cxd 
and the system will be ordered at sufficiently low T and 
one expects a phase transition at some Tc > between 
two distinct phases. The standard method to implement 
these ideas is to simply find the energies of the system 
subject to periodic and antiperiodic BC, despite the fact 
that neither BC is compatible with the true GS and that 
the energy difference is not a defect energy but is the 
energy difference between two randomly chosen defects. 
Nevertheless, < AE{L) >=< \Eap{L) - Ep{L)\ > is as- 
sumed to scale as and the best fit to this is taken as 
yielding the exponent 9. 

In this paper, we take the point of view that, since ran- 
domness restricts the accessible system sizes L to rather 
small values, for the DWRG to be successful it is im- 
portant to minimize uncontrolled and poorly understood 
effects which may introduce large corrections to scaling. 
Also when fitting a continuous function to less than ten 
data points there is scope for obtaining erroneous results. 
To attempt to minimize such effects, we first transform 
the problem described by eq. (|l|) to a Coulomb gas repre- 
sentation which eliminates the spin wave contribution to 
the energies and also enables us to estimate numerically 
the true GS energy by optimizing the BC. In the phase 
representation this would involve applying phase differ- 
ences across all the bonds beween opposite faces of the 
system on a hypertorus in d dimensions and minimizing 
the energy with respect to L'' phases and dL'^^^ phase 
differences on the bonds connecting opposite faces. In 
the Coulomb gas representation, one needs only d extra 
parameters, one for each direction, and minimize with 
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respect to {d + L"^) parameters. To introduce a domain 
wall, we have to change only one of the d parameters 
from its value in the GS and minimize the energy with 
respect to L'^ bulk variables keeping the BC fixed to find 
the defect energy Ed. We then obtain the true domain 
wall energy as AE{L) = Ed{L) — Eo{L) > for a partic- 
ular realization of disorder. This procedure is repeated 
for several realizations of disorder and system sizes L to 
obtain < AE{L) > which is fitted to eq.(||). 

Some support for this point of view comes from a con- 
jecture made by Ney-Nifle and Hilhorst [Q in the con- 
text of the 2d XY spin glass, which is a special case of 
the model of eq.(|l|) when the bond variables Aij = 0, tt 
each with probability 1/2, based on some analytic work 
by Ney-Nifle, Hilhorst and Moore on the XY spin 
glass on a ladder. The latter authors solve the problem 
analytically and, by imposing reflective BC which forces 
a chiral domain wall into the system, find that both the 
spin and chiral domain wall energies scale with the same 
exponent 6s = Oc = — 1.7972 The same conclusion 
is reached in ref. but with a slightly different value 
of 9. In ref. |0, supported by heuristic but plausible 
arguments, it is conjectured that, for any spin glass sys- 
tem below the lower critical dimension di, 9s = 6c < 0. 
This result, although lacking rigorous proof, provides an 
important bench test for a numerical simulation which 
seems to have been either ignored or overlooked in some 
recent studies eg. [0 where it is found that 9s < 9c in 
d = 2 < di. On the other hand, we have recently stud- 
ied the 2d XY spin glass using the philosophy outlined 
above and find that 9s = 9c ^ —0.35 to the accuracy 
of our simulations p8[ |. This supports the point of view 
taken in this paper that defect energies should be mea- 
sured by changing the BC in a controlled way from those 
BC consistent with the true ground state of a system 
with a particular realization of randomness. If this is not 
done, as in the standard method of comparing energies 
with periodic and antiperiodic BC, one is not measuring 
a specific defect energy and the < AE{L) > obtained 
this way either does not scale as or 9 is the stiffness 
exponent of a quantity whose meaning is unclear. In any 
event, if the latter scenario holds, there is no reason to 
expect any relation between 9s and 9c. Since the XY 
spin glass is a special case of the gauge glass, we expect 
that these considerations also hold for the latter system. 

The same cannot be said for the standard RT method 
where there are no checks, just a hope, that true minima 
have been achieved. To eliminate the effect of smooth vari- 
ations in the phases 9i, we transform the Hamiltonian 
from the phase representation of cq.(|l|) to a Coulomb 
gas representation which is more convenient for numeri- 
cal work. The function V{(p) is taken to be a piecewise 
parabolic potential which is equivalent to a Villain | [T9| ] 
potential at T = 
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(3) 



By a duality transformation 21 1, in 2d the 

Coulomb gas Hamiltonian becomes 



2 7rV5](gr-/r)G(r-r')(gr'-/r') 

r,r' 

2Z^ 51 



(4) 



a—x,y 



where 



(Jx = —211 



(jy — — 27r 



L{qxi - fxi) + ^{qv- fr)y 

r 

L{qyl - fyl) - ^{Qr - fr)x 



N ^ 4 — 2 cos kx — 2 cos fc„ 



(5) 



r = {x,y) represents the coordinates of the sites of the 
dual lattice and G(r) is the lattice Green's function. The 
topological charge, Qr, is the circulation of the phase 
round the plaquette at r, — X^Dr'^y/^"' ^'^'^ 
any integer subject to the neutrality condition q^. = 0. 
The frustration, /r, is the circulation of Aij/27r around 
the plaquette at r. The quantity jx\ is the circulation 
of Aij /2tt round the whole torus on horizontal bonds of 
the plaquettes at y = 1 and qxi is the coresponding cir- 
culation of the phase, fyi and qyi are defined similarly. 
Periodicity in the phases 9i restricts qxi,qyi to integers. 
In 3(i, the charge Hamiltonian becomes 



H = 2 7rVE(qr - fr) • (qr' - fr')G(r - r') 

r.r' 



(6) 



where 



ax = TtL-' J2{-<l'r - -fr) + yiQr " fr)} + Q. 
r 

Qx = -nY^i-zfySy^i+yrj^^ + 2TiL{qxi - fxi) 



G(r) = 1 y _ 



N ^ 6 — 2 cos kx — 2 cos ky — 2 cos kz 

kyiO X y z 



(7) 



with similar expressions for cr^ , CTz , Qy , Qz which are ob- 
tained by cyclic permutations of (x, y, z). Here, as in 2d, 
r = (a;, y, z) are the coordinates of the dual lattice sites 
and G(r) is the Green's function. The charge qr and the 
frustration fr associated with the dual lattice at r be- 
come vector quantities with components {q^,q^,q^) and 
{fr ^ fr ^ fr)- Here and /" are the circulations of 9i 
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and Aij about the plaquette at r normal to the direction 
a. These vector charges satisfy V • qr = at each site r 
and a neutrahty condition as in 2d. 

To find the true GS energy of the finite system of L'^ 
sites on a hypertorus, we minimize the energy given by 
eq.(Q) with respect to the variables qr and fxi, fyi in 2d, 
and by eq.(^ in 3d. The lowest energy of the system 
is 27r periodic in the twist = 27r/^i, with a mini- 
mum at some A° which depends on the particular real- 
ization of disorder. The lowest energy with ^ Aj^ 
includes the excitation energy due to the twist A^ — A" . 
Adding a twist is equivalent to the gauge transformation 
Aij Aij + A^/L on all bonds in the direction /i. 

We compute the domain wall energy AE{L) in two 
different ways. The first is with the usual periodic and 
antiperiodic BC, which correspond to twists A ^ A° 
and A -I- ttx where A is determined by the particular 
sample which we call a random twist (RT) measurement. 
The second is by computing the true domain wall energy 
AE^'^iL) = El{A° + ttx) - El{A°) where £;l(A°) is 
obtained by minimizing the energy with respect to the 
bulk charges and the twists fia- The domain wall 
is induced by fi^ -1-1/2 which corresponds to 

A^ A|^ -I- TT which we call a best twist (BT) mea- 
surement. Note that, since El{AF') is the absolute en- 
ergy minimum, AE^'^{L) > and the boundary con- 
tributions to the total energy proportional to cr„ > 
in eqs.(^j^) must vanish when the best twist A° is ap- 
plied. These may be used as checks on the numerical 
algorithm. We use simulated annealing [^,^ to esti- 
mate the lowest energies, which is considerably more ef- 
ficient than repeated quenches to T = 0. To minimize 
errors in the computation of domain wall energies, we 
used two random number sequences to determine the an- 
nealing schedule and demanded that both sequences yield 
the same GS energy, except for our largest system in Sd 
with L = 7 because of time constraints. The main source 
of error is statistical due to the averaging over disorder, 
but the ratio AE{L)/ El{A°) decreases rapidly with in- 
creasing L, especially in 3d where this ratio is about 10^^ 
for L = 7, which means that £'l(A) must be obtained 
almost exactly. This puts severe constraints on the ac- 
cessible system sizes L. 

In 2d, we compute the domain wall energy for sizes 
2 < i < 10, taking averages over about 10'^ realiza- 
tions of disorder for each L. The results of RT and 
BT measurements are shown in fig.(l). The RT mea- 
surement gives O^'^ = —0.45 ± 0.015 which agrees with 
other groups, all of whom have used this measurement 
[^-|| . This agreement is not surprising as the only differ- 
ence between these is in the computing power and in the 
algorithms used. All are measuring the same quantity 
A£'^-^(i). The true defect energy obtained from the BT 
measurement gives 9^'^ = —0.36 ±0.010 which is consid- 
erably larger. In 3d, the system sizes are 2 < L < 7 with 
disorder averaging over 10"^ realizations for L < 5, 300 



for L = 6 and 60 for L = 7. The error in AE{L = 7) is 
very large, but this point was included to check that it 
is consistent with the behavior deduced from the smaller 
systems. The results are shown in fig. (2). For the RT 
measurement, there seems to be a crossover around L = 5 
from a very small value of 9^'^ to a larger positive value as 
is also seen in |^ , but our sizes do not allow any estimate 
of 9^^. On the other hand, the BT measurement for 
these sizes is consistent with a power law scaling with a 
stiffness exponent 9^'^ = +0.31 ±0.010 using sizes L < 6, 
which is strong evidence in favor of a superconducting 
glass phase at finite T and of d; < 3. This is also consis- 
tent with finite T Monte Carlo results on the 3d gauge 
glass P,^,^ which indicate Tc ~ 0{J), which is difficult 
to reconcile with the very small value of 9^^ which, if it 
were the stiffness exponent for the 3d gauge glass, would 
imply d; w 3 and a small value of Td J . If we include 
the L = 7 point, the best fit gives 9^^ = +0.30 ± 0.015. 
Note that the errors quoted here in Q^t,rt ^Qj^g from 
a naive least squares fit to the data, and should not be 
taken too seriously. The L = 7 data is suspect because, 
in 10^ CPU hours on a Cray J90, 4 samples of a batch 
of 64 violated the BT condition AE^'^{L) > implying 
insufficient annealing to reach the true energy minima. 
What data we have is entirely consistent with the scaling 
form of eq.(H) with 9 w +0.3 with no sign of any devia- 
tion from this. 

We also studied the effects of screening on the domain 
wall energy using the BT measurement in 3d. Screening 
of the interaction of charges is implemented by adding a 
term A^^, where A is the screening length, to the denom- 
inators of the Green's function of eq. (|^) |^ . The results 
are also shown in fig. (2). We averaged over 10"^ samples 
for L = 2, 3, 4 and 250 for L = 5. For the shorter screen- 
ing lengths, screening is clearly a relevant perturbation 
and destroys the ordered phase while for longer screen- 
ing lengths AE{L) seems to scale the same way as the 
unscreened case but we expect there is a crossover to a 
negative stiffness exponent at length scales L which are 
inaccessible with our computing power. Our results are 
consistent with those of Bokil and Young who stud- 
ied the screening question using the RT measurement. 

The major result of this study is that one can, in prin- 
ciple, find the exact GS energy of a random XY system 
by a suitable choice of boundary conditions which are 
consistent with the unknown GS. This is implemented 
for the gauge glass model in both 2d and 3d and it is 
also argued that a domain wall is created from the GS 
by an appropriate change of the BC. In the Coulomb gas 
representation of the system on a hypertorus, the BC are 
parametrized by d numbers which are the circulations of 
Aij round loops enclosing the hypertorus. Once these are 
known, a domain wall is induced by changing one of these 
by TT. The size dependence of the domain wall energy is 
computed and is found to scale rather accurately as L^. 
The values of the stiffness exponent 9 are much larger 
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than previous estimates and is in accord with < 3. 
This also reconciles the finite T Monte Carlo results with 
one's physical intuition. We stress that the disagreement 
between the stiffness exponent 9^^ in this work and all 
previous estimates is because these measure O^'^ which 
is a quantity whose meaning is unclear and is more likely 
to suffer from large corrections to scaling, especially for 
the small L values which can be simulated (see fig. (2) 
and ref. ||]). We conjecture that they would coincide if 
much larger values of L could be reached. The stiffness 
exponent 9^^ is larger than the 9^^'-^ « +0.2 for the 3d 
Ising spin glass p6| but we see no contradiction here as 
continuous variables can adjust to frustration more eas- 
ily than Ising spins and one expects gauge glass order to 
resist distortions better than Ising spin glass order. 

Computations were performed at the Theoretical 
Physics Computing Facility at Brown University. JMK 
thanks A. Vallat and B. Grossman for countless conver- 
sations on gauge glasses when many of the ideas in this 
paper were proposed. 



[1] 

[2] 



[4; 

[5] 
[6] 

[7 

[8] 
[9] 

[lo; 
[11 

[12 

[13 
[14 

[15 

[16 

[17 

[18 



M. P. A. Fisher, Phys. Rev. Lett. 62, 1415 (1989) 

G. Blatter, M. V. Feigel'man, V. J. Geshkenbein, A. I. 

Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 

(1994) 

D. A. Huse and H. S. Seung, Phys. Rev. B 42, 1059 
(1990) 

J. D. Reger, T. A. Tokuyasu, A. P. Young, and M. P. A. 

Fisher, Phys. Rev. B 44, 7147 (1991) 

M. J. P. Gingras, Phys. Rev. B 45, 7547 (1992) 

M. P. A. Fisher, T. A. Tokuyasu, and A. P. Young, Phys. 

Rev. Lett. 66, 2931 (1991) 

J. M. Kosterlitz and M. V. Simkin, Phys. Rev. Lett. 79, 
1098 (1997) 



J. Maucourt and D.R. Grempel, |cond-mat/9802242 
C. Dekker and P. J. M. Woltgens, R. H. Koch, B. W. 
Hussey, and A. Gupta, Phys. Rev. Lett. 69, 2717 (1992) 
R. H. Koch, V. Foglietti, W. J. Gallagher, G. Koren, A. 
Gupta, and M. P. A. Fisher, Phys. Rev. Lett. 63, 1511 
(1989) 

C. Dekker, W. Eidelloth, and R. H. Koch, Phys. Rev. 
Lett. 68, 3347 (1992) 

J. R. Banavar and M. Cieplak, Phys. Rev. Lett. 48, 832 
(1982) 

W. L. McMillan, Phys. Rev. B29, 4026 (1983) 

M. Ney-Nifle and H. J. Hilhorst, Phys. Rev. B51, 8357 

(1995) 

M. Ney-Nifle, H.J. Hilhorst and M.A. Moore, Phys. Rev. 
B48, 10254 (1993) 

M.J. Thill, M. Ney-Nifle and H.J. Hilhorst, J. Phys. A28, 
4825 (1995) 

J. Maucourt and D. R. Grempel, Phys. Rev. Lett. 80, 
774 (1998) 

J.M. Kosterlitz and N. Akino, unpublished 



[19] J. Villain, J. Phys. (Paris) 36, 581 (1975) 

[20] J. V. Jose, L. P. Kadanoff, S. Kirkpatrick, and D. R. 

Nelson, Phys. Rev. 16, 1217 (1977) 
[21] A. Vallat and H. Beck, Phys. Rev. B50, 4015 (1994) 
[22] C. Wengel and A.P. Young, Phys. Rev. B56, 5918 (1997) 
[23] S. Kirkpatrick, CD. Gellat and M.P. Vecchi, Science 

220, 671 (1983) 
[24] M.V. Simkin, Phys. Rev. B55, 11405 (1997) 
[25] H. S. BokU and A. P. Young, Phys. Rev. Lett. 74, 3021 

(1995) 

[26] A.J. Bray and M.A. Moore, J. Phys. C17, L463 (1984); 
W.L. McMillan, Phys. Rev. B31, 340 (1985) 



S. 0.5 - 



0.0 





O^BT( 


nscreened) 






nscreened) 









0.5 



1.0 



1.5 
ln{L) 



2.0 



2.5 



FIG. 1. Size L dependence of domain wall energy in 
2d. Both RT and BT measurements are shown. Solid 
lines are power law flts. Error bars are not shown if 
smaller than symbol size. 
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FIG. 2. L dependence of domain wall energy in Ttd. 
Bottom curve is RT measurement for unscreened inter- 
action. All others are BT measurements. Topmost curve 
is unscreened case L = 2 — 7. Other curves are screened 
interactions with A decreasing from top to bottom. Solid 
lines are power law fits. Dotted lines are guides for the 
eye. 
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